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The properties of two-component Fermi gases with zero-range interactions are universal. We use 
an explicitly correlated Gaussian basis set expansion approach to investigate small equal-mass two- 
component Fermi gases under spherically symmetric external harmonic confinement. At unitarity, 
we determine the ground state energy for systems with up to ten particles interacting through finite- 
range two-body potentials for both even and odd number of particles. We extrapolate the energies 
to the zero-range limit using a novel scheme that removes the linear and, in some cases, also the 
quadratic dependence of the ground state energies on the two-body range. Our extrapolated zero- 
range energies are compared with results from the literature. We also calculate the two-body Tan 
contact and structural properties. 
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I. INTRODUCTION 

The properties of two-component Fermi gases interact¬ 
ing through two-body zero-range potentials with s-wave 
scattering length a s are universal [1, 2]. At unitarity, 
i.e., for infinitely large a s , the two-body interaction does 
not define a meaningful length scale and the strongly- 
interacting Fermi gas is characterized by the same num¬ 
ber of length scales as the non-interacting Fermi gas. Ap¬ 
proximate realizations of the unitary Fermi gas include 
dilute neutron matter in the crusts of neutron stars [3] 
and ultracold atom gases such as 6 Li [4] and 40 K [5]. The 
properties of homogeneous and inhomogeneous unitary 
Fermi gases have attracted a great deal of experimen¬ 
tal and theoretical attention. For spherically symmetric 
external confinement, the harmonic oscillator length aho 
defines the only length scale of the system. It is hence 
interesting to determine how the properties of trapped 
unitary Fermi gases vary with the number of particles. 

Harmonically trapped Fermi gases at unitarity have 
been treated by quantum Monte Carlo methods [6-13], 
density functional theory (DFT) [11, 14-17], and basis set 
expansion approaches. The accuracy of the fixed-node 
diffusion Monte Carlo (FN-DMC) energies [7, 8, 11, 12] 
depends on the quality of the many-body nodal surface. 
The resulting energies provide upper bounds to the exact 
ground state energies and the zero-range limit is reached 
through extrapolation [11, 12]. Auxiliary-field quantum 
Monte Carlo (AFMC) methods, on the other hand, work 
on a finite lattice and extrapolation to the infinite lattice 
limit is required to obtain fully converged result [12]. The 
quality of DFT calculations depends critically on the un¬ 
derlying functional. Since the functional is typically ob¬ 
tained by matching to data for the homogeneous system, 
the analysis of results for the trapped system can provide 
insights into gradient corrections and other finite size fea¬ 
tures [15-17]. Trapped unitary Fermi gases with up to 
six particles have been calculated by explicitly correlated 
Gaussian (ECG) basis set expansion approaches [7, 8, 18- 
22] with better than about 1% accuracy. Application of 
the ECG method to systems with more than six particles 


has been challenging due to the rapid increase of the num¬ 
ber of permutations and the larger number of degrees of 
freedom. Recently, Ref. [23] treated the (TVi,TV 2 ) = (4,4) 
system at unitarity using a basis set that accounts for the 
most important but not all correlations. 

Here, we present results for small trapped unitary 
Fermi gases with TV < 10 particles, where TV = + TV 2 

and N\ — TV 2 = 0 or 1. Our extrapolated zero-range en¬ 
ergy of the (4,4) system is 0.9% lower than that reported 
in Ref. [23]. A new aspect of our work is that we de¬ 
veloped an improved scheme for extrapolating the finite- 
range energies to the zero-range limit. This new scheme 
eliminates the linear and, in some cases, the quadratic 
dependence of the ground state energies on the two-body 
range. The scheme provides a consistency check on the 
range-dependence of our energies and reduces the errors 
that result from the extrapolation to the zero-range limit. 
Our results suggest that the developed range correction 
scheme allows one to obtain a reliable approximation to 
the zero-range energy from a single finite-range calcu¬ 
lation. The scheme can be applied to other numerical 
calculations that work with finite-range interactions. We 
use our range correction scheme to determine the zero- 
range energies and the Tan contact for two-component 
Fermi gases with TV < 10 at unitarity. In addition, we 
present selected structural properties. 

The remainder of this paper is organized as follows. 
Section II discusses the theoretical framework and our 
extrapolation scheme to the zero-range limit. Section III 
presents our results for systems with up to ten particles 
and compares, where available, with results from the lit¬ 
erature. Lastly, Sec. IV concludes. 


II. THEORETICAL FRAMEWORK 

We consider equal-mass two-component Fermi gases 
with N-\ spin-up and TV 2 spin-down atoms (TV = Aq + 
TV 2 and TVi — TV 2 = 0 or 1) under external spherically 
symmetric harmonic confinement with angular trapping 
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frequency to. The system Hamiltonian H(r 0 ) reads 

N n 2 

H (r 0 ) = E - -Sn) 

2—1 

Ni N 

+E E V 2 b(rij,r 0 ), (1) 

2—1 j = Ni~\~l 

where m denotes the atom mass, r, denotes the position 
vector of the ith particle with respect to the trap center, 
and 

N 1 

Vt r (ri,..-,fjv) = E 2 TOw2 ^ ( 2 ) 

i—1 

is the trapping potential. V 2 b is the interspecies two-body 
interaction potential that depends on the interparticle 
distance rq, rq = r) — fj\. In our work, it is modeled 
by a finite-range Gaussian potential with range rn and 
depth E 0 (Vb < 0), 

V 2 b(r, r 0 ) = V 0 exp ^ . (3) 

For a fixed r 0 , Vq is adjusted such that V 2 b(r, r 0 ) has an 
infinitely large s-wave scattering length a s and supports 
one zero-energy two-body bound state in free space. The 
ranges r 0 considered depend on the size of the system 
and vary from O.Olaho to 0.12oh o , wh ere a^ 0 d enotes the 
harmonic oscillator length [aho = ^/S/(mw)]. For the 
Gaussian potential with one zero-energy bound state, the 
effective range r e s is approximately equal to 2.032ro. 

To numerically solve the Schrodinger equation for the 
Hamiltonian given in Eq. (1), we separate off the cen¬ 
ter of mass degrees of freedom and expand the eigen¬ 
states of the relative Hamiltonian in terms of ECG ba¬ 
sis functions, which depend on a set of non-linear vari¬ 
ational parameters that are optimized through energy 
minimization (see below) [18, 24]. The unsymmetrized 
basis functions for states with L 7r = 0 + and 1“ sym¬ 
metry (L denotes the relative orbital angular momen¬ 
tum and 7T the relative parity) read exp(— ^x T Ax) and 
y\o{u T x) exp(— \x T Ax), respectively, where A is a sym¬ 
metric and positive definite ( N — 1) x (N — 1) parameter 
matrix, u = (iq, u 2 , un-i) t is a IV - 1 dimensional 
vector, and J4o is a solid spherical harmonic function [24], 
x = (xi, x 2 , ■■■,xn-i) t collectively denotes a set of N — 1 
Jacobi vectors. The ground state of even N systems has 
0+ symmetry and that of odd N systems has 1“ symme¬ 
try. A key advantage of these basis functions is that the 
corresponding overlap and Hamiltonian matrix elements 
can be calculated analytically [24]. 

The fermionic exchange symmetry is ensured by acting 
with the permutation operator A on the unsymmetrized 
basis functions. The number of permutations N p in¬ 
creases factorially with the number of identical fermions. 
For the (5, 5) system, e.g., A contains (5!) 2 = 14,400 two- 
particle exchange operations with alternating plus and 


minus signs. The evaluation of each overlap and Hamil¬ 
tonian matrix element involves a sum over N p terms that 
are highly oscillatory. In a standard 16 digit floating 
point implementation, numerical challenges arise from 
the near-cancellation of the positive and negative terms 
for systems with N > 8. The near-cancellation of these 
terms of alternating signs can be interpreted as a rel¬ 
ative of the fermion sign problem known from Monte 
Carlo simulations [25, 26]. To ensure that the matrix 
elements for the largest systems considered are accurate 
to at least ten significant digits, we implemented our C 
codes using extended precision. The eigenenergies and 
expansion coefficients are obtained by solving a gener¬ 
alized eigenvalue problem that involves the Hamiltonian 
matrix and the overlap matrix. The numerical error of 
the resulting eigenenergies is several orders of magnitude 
smaller than the errors that arise from the use of a finite 
basis set and the extrapolation to the zero-range limit. 
In Sec. Ill, we report the total ground state energy E(ro) 
of the Hamiltonian H(ro), i.e., we add the center of mass 
energy of 3Eh 0 /2 to the relative energy obtained by the 
ECG approach. Here, E-o 0 denotes the harmonic oscilla¬ 
tor energy (Eho = fuS). 



FIG. 1. (Color online) Illustration of the Jacobi coordinates 
employed in our work for the (Ni,N 2 ) = (5,4) system. The 
dark vectors show the Jacobi vectors X\. x 2 , ..., xs. The spin- 
up and spin-down fermions are represented by light vertical 
up and down arrows. 

We use a senri-stochastic variational approach to 
choose and optimize the variational parameters contained 
in A and u [24]. Our Jacobi coordinates are chosen such 
that the first N 2 Jacobi vectors correspond to distance 
vectors between unlike particle pairs. The next -/V 2 /2 Ja¬ 
cobi vectors correspond to the distance vectors between 
the center of mass of the first pair and the second pair, the 
distance vector between the center of mass of the third 
pair and the fourth pair, and so on. The remaining Ja¬ 
cobi vectors connect the larger sub-units and, for odd N, 
the iVth particle (see Fig. 1 for an illustration for N = 9). 
For this choice of Jacobi coordinates, the first N 2 diag¬ 
onal elements of the A matrix represent correlations be¬ 
tween unlike particles. We expect that these interspecies 
distances are, on average, smaller and more strongly cor¬ 
related than those between like particles. This motivates 
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us to choose the first _ZV 2 diagonal elements of the A ma¬ 
trix for each basis function from a preset non-linear grid. 
The other diagonal elements are chosen stochastically 
from preset parameter windows. As in Ref. [23], we start 
with basis functions that are diagonal in A. These basis 
functions account for the most important correlations. 
Once a certain basis set size is reached (for the larger 
systems around Nb = 500), we reoptimize the variational 
parameters contained in the diagonal of A , and for odd 
N in u, and allow for off-diagonal A matrix elements. 
The off-diagonal matrix element at position (i,j) of the 
A matrix is chosen from 0 to the geometric mean of the 
zth and j th diagonal elements. We find that this choice 
results with high probability in positive-definite A matri¬ 
ces. The positive-definiteness of the A matrix is tested 
through diagonalization. We refer to the reoptimization 
of all variational parameters contained in A and, if appli¬ 
cable, u of all Nb basis functions as a reoptimization cy¬ 
cle. The re-optimization cycle is repeated until the lowest 
energy changes by less than a preset value. After that, 
we extend the basis set by a hundred to several hundred 
basis functions and reoptimize the variational parame¬ 
ters of the enlarged basis set using a variable number 
of reoptimization cycles. This process is repeated a few 
times. At the end, the basis is enlarged to around 2000 
basis functions without additional reoptimization of the 
non-linear variational parameters. The basis set errors 
reported in Tables I and II of Sec. Ill and Tables I-VI 
of the Supplemental Material [27] are estimated by ana¬ 
lyzing the energy decrease that results from the basis set 
enlargements and the reoptimization cycles. 

To reach the universal regime where aho defines the 
only length scale in the system, we need to extrapolate 
the numerically calculated finite-range energies to the 
zero-range limit. In previous ECG works [18, 19, 21 
23], this was done by fitting the finite-range energies by 
a linear or quadratic function. We refer to this tra¬ 
ditional extrapolation scheme as the zeroth-order ex¬ 
trapolation scheme. The difference between the finite- 
range energies and the extrapolated zero-range energies 
is, typically, at the order of a few percent and can intro¬ 
duce a non-negligible extrapolation error. Moreover, for 
larger systems, it is computationally expensive, maybe 
even prohibitively expensive, to obtain energies at very 
small ranges. It should also be noted that the extrapo¬ 
lated zero-range energies do not provide variational up¬ 
per bounds even though the finite-range ECG energies 
do. It is thus desirable to remove the linear and, ideally, 
quadratic range dependence. Motivated by the general¬ 
ized virial theorem 


E{ 0) = 2V tr (0) (4) 

[ktr(0) denotes the expectation value of Vt r (fi, fjv) for 
r 0 —> 0] at unitarity, Werner [28] proposed to remove the 
linear range-dependence of the ground state energy by 
combining it with the expectation value Vt r (ro) of the 
trapping potential Vt r (ri, ...,rjv) calculated for the same 


ro, 

E{ 0) = 3 E(r 0 ) - 4Vtr(r 0 ) + 0(rg). (5) 

While Eq. (5) removes the leading-order range depen¬ 
dence, it is associated with errorbars that come from the 
basis set errors of E(r 0 ) and Vtr(^o)- In our ECG method, 
the basis set is optimized by minimizing the ground state 
energy. Not surprisingly, we find that the convergence of 
the expectation value of the trapping potential is not as 
good as that of the energy. This motivates us to propose 
an alternative scheme that can be carried out to higher 
orders. 

The ground state energy E(ro) of the IV-particle sys¬ 
tem is a smooth function of the two-body interaction 
range rg. The n max th order Taylor series of E(fo) around 
r 0 is 


Umax -i 

E(f 0 )= VB (n) (ro)w(fo-r 0 r 

L ' 77 , 


n—0 


+O[(r 0 - r 0 ) 


71 max “1“ 1] 


where 


E {n \r 0 ) = 


d n E(f 0 ) 

<9fg 


( 6 ) 


( 7 ) 


is the nth order derivative of the ground state energy with 
respect to the range evaluated at r 0 . E^\ro) is simply 
the ground state energy E(ro) of H(r 0 ). E^’{r o) can be 
obtained through the Hellmann-Feynman theorem [29], 


e (1) M 


/dH(f 0 ) 
\ df 0 


( 8 ) 


which is exact in the limit that the basis set is complete. 
The matrix elements needed to evaluate E^\ro) reduce 
to compact analytical expressions. E^ 2 \ro) can be ob¬ 
tained by the finite difference method, i.e., by evaluating 
E^(ro) at two nearby r 0 . 

Our goal is to obtain the zero-range energy E(0). Set¬ 
ting ro in Eq. (6) to 0, we obtain 

E( 0) = E ZRA ,„ max (ro) + 0(rg max+1 ), (9) 

where 

77 max -i 

£zRA,n max (ro) = £ E(n) ( r °)~i (-ro)" (10) 

' nl 

n=0 


is the JZmaxth order approximation to the zero-range en¬ 
ergy. Equations (9) and (10) establish a relation be¬ 
tween the zero-range energy E(0) and the finite-range 
energy E(r 0 ) and its derivatives with respect to the two- 
body range. EzRA,o( r o) is simply the finite-range en¬ 
ergy E(ro) with linear leading-order range-dependence. 
The leading-order range-dependence of -EzRA,i(fo) and 
Ezra, 2 (^ 0 ) is quadratic and cubic, respectively. Cru¬ 
cial is that ^zra,i(?*o) and Ezra, 2 (^ 0 ) are obtained 
at finite ro without extrapolation. They provide bet¬ 
ter approximations to the zero-range energy E{ 0) than 
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-EzRA,o( r o)- We refer to the extrapolations of E Z R A .i(r 0 ) 
and -Ezra,2(^0) to the zero-range limit as the first- and 
second-order extrapolation schemes. For a complete ba¬ 
sis, Ezra,1(^0) coincides with the quantity 3 E(ro) — 
44tr( r o)i he., formally Eq. (9) with n max = 1 is equiv¬ 
alent to Eq. (5). It turns out, however, that our ECG 
implementation provides a more accurate estimate for 
E^(?'o) than for Vtr(^o)- In Sec. Ill, we independently 
fit Ezra,0O0), Ezra.i^o), and E Z R A , 2 (r 0 ) and compare 
the resulting zero-range energies. Appendix A shows 
that the functional forms of Ezra,o(a>)> Ezra, i(ro), and 
Ezra,2(^0) are correlated and presents the results of 
a single combined fit of E Z R Ai o(r 0 ), E Z R A ,i(r 0 ), and 
Ezra,2(70)- The resulting zero-range energy is found to 
be consistent with the zero-range energies obtained from 
the independent fits. 

The ECG calculations become numerically more chal¬ 
lenging with decreasing two-body interaction range ro- 
The challenges arise from the need to resolve length scales 
of different orders of magnitude. In previous ECG cal¬ 
culations [21, 22], much effort was put on solving the 
Schrodinger equation for systems with small r 0 . In this 
work, we show that a reliable approximation to the zero- 
range energy can be obtained by calculating Ezra, 2 (^ 0 ) 
at a range ro « O.laho- We demonstrate in Sec. Ill that 
EzRA.i( r o) and Ezr A , 2 (To) play important roles in ob¬ 
taining the zero-range energy E( 0) and its errorbar. 

To calculate the two-body Tan contact C(ro) at uni- 
tarity, we use the adiabatic energy relation [30], 


C(r 0 ) 


47tto <9E(ro) 

h 2 d^—aj 1 ) 



( 11 ) 


To obtain the two-body contact for ro = 0, we use the 
zeroth- and first-order zero-range extrapolation schemes, 
i.e., we extrapolate Czra,o(?'o) and 

ro (12) 

ro~* r 0 

to the zero-range limit. 

The contact can alternatively be calculated through 
the pair relation 

C(r 0 ) = Ni x N 2 x lim (47r) 2 Pi 2 (r)r, (13) 

r—>0,r^$>ro 


Gzra.iM = C'(ro) - PP 


where Pi 2 (r) denotes the pair distribution func¬ 
tion. The quantity r 2 Ei2(r) with normalization 
4 tt f 0 °°P 12 (ry dr = 1 tells one the likelihood of find¬ 
ing two unlike particles at distance r from each other. 
The behavior of 47rPi2(r)r 2 around r ss ro depends on 
the details of the two-body interaction potential. Specif¬ 
ically, for finite-range potentials the quantity 47rP 12 (r)r 2 
goes to zero as r — > 0. For the zero-range potential, in 
contrast, 47rPi2(r)r 2 remains finite as r — > 0. Thus, to 
extract the finite-range contact via the pair relation, we 
consider the region where r r 0 but r <C Oho- In prac¬ 
tice, the condition r>ro translates to r > 2r 0 

We also consider the spherically symmetric radial den¬ 
sity Pj(r ) of species j. j = 1 and 2. For even N, 


we have Pi(r) = P 2 (r). The quantity Pj(r) tells one 
the likelihood of finding a particle at distance r from 
the trap center. The normalization is chosen such that 
47r /o°° Pj(r)r 2 dr = 1. 


III. RESULTS 



FIG. 2. (Color online) Ground state energy of the (3,2) 
system at unitarity as a function of r 0 . Circles, squares, 
and diamonds show the energies EzRA,o(ro), EzRA,i(ro), and 
Ezra,2(^0), respectively, reported in the last three columns 
of Table I. The uncertainty of AE(ro) is not accounted for by 
the errorbars. Solid, dashed, and dotted lines show polyno¬ 
mial fits to E Z RA,o(ro), E Z ra,i(?~o), and E Z RA,2(r 0 ). 

This section discusses the energies and other observ¬ 
ables of the systems with N < 10 obtained by the ECG 
method. We use the (3,2) system to explain our new 
range correction scheme. Column 2 of Table I shows 
the finite-range energies E(ro) obtained by the ECG ap¬ 
proach for various two-body interaction ranges ro- The 
reported energies are obtained for the largest basis set 
considered. They provide variational upper bounds for 
the finite-range Hamiltonian with Gaussian interaction. 
Column 3 reports the estimated basis set error AE(ro). 
For all ro considered, the basis set error is less than 
0.02%. Columns 4 and 5 show the quantities E^ (r 0 ) and 
E^ 2 ^(ro), respectively. While EW(r 0 ) increases slightly 
with decreasing range r 0 , this increase is smaller than 
the decrease of r 0 , implying that the range correction 
E- 1 - 1 (ro)r'o decreases with decreasing ro- The magnitude 
of E^ 2 \ro) decreases with decreasing ro- Note that we 
are not able to estimate E^ 2 \ro) reliably for small ro 
(the errorbars are larger than the quantity itself). Yet, 
the errorbars of E^ 2 \ro) allow us to estimate the max¬ 
imal correction proportional to r^ for each ro, thereby 
providing us with another means to estimate errorbars. 
Columns 6-8 of Table I show E Z R A j(r 0 ) with j = 0, 
1, and 2. These values are obtained by subtracting the 
basis set error AE(r 0 ). The leading-order range depen¬ 
dence of EzRA,o( r o) is linear and we perform a fit of the 
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TABLE I. Ground state energy of the (3,2) system at unitarity. Column 2 shows the finite-range energy for the largest basis 
set considered. The estimated basis set error AE(ro) is reported in column 3. Columns 4 and 5 report the quantities E^(ro) 
and E^ 2 \r o); errorbars are given in parenthesis. The energy derivatives are calculated for the largest basis set considered. 
Columns 6-8 report the energies Ezra,o()"o), EzRA,i(ro), and Ezra, 2 (^ 0 ). These energies account for the estimated basis set 
extrapolation error, i.e., E(r 0 ) — AE(ro) is being used to calculate EzRA,j(ro) for j = 0, 1, and 2. The errorbars of EzRA,i(ro) 
and Ezra, 2 (^ 0 ) account for the uncertainties of E^(ro) and E (2 \r 0 ) but do not account for the uncertainty of AE(ro). The 
last row reports the extrapolation of EzRA,j(ro) to the zero-range limit. 


.roi. E ( r o) 

AE(r 0 ) 

E (1) (ro )^t 


■Ezra,0(^0) 

Ezra, 1(^0) 

E hn 

Ezra, 2(^0) 
Eho 

0.07 7.5449 

0.0061 

1.148(7) 

-5.13(12) 

7.5448 

7.4644(5) 

7.4518(8) 

0.06 7.5332 

0.0001 

1 . 201 ( 8 ) 

-4.61(25) 

7.5331 

7.4610(5) 

7.4527(9) 

0.05 7.5211 

0.0004 

1.240(12) 

-4.11(28) 

7.5207 

7.4587(6) 

7.4536(10) 

0.04 7.5084 

0.0004 

1.284(15) 

-2.85(33) 

7.5080 

7.4566(6) 

7.4543(9) 

0.03 7.4954 

0.0006 

1.297(26) 

-0.32(46) 

7.4948 

7.4559(8) 

7.4557(10) 

0.02 7.4825 

0.0009 

1.290(120) 

0.24(96) 

7.4816 

7.4558(24) 7.4559(26) 

0 




7.4557 

7.4550 

7.4563 


form co+Ciro + C 2 rg-)-C 3 rQ. The extrapolated zero-range 
energy is reported in the last row of column 6. The lead¬ 
ing order range dependence of Ezra, 1 (no) is quadratic 
and we perform a fit of the form Co + C 2 rg + C 3 Tq + c^Tq, 
weighted by the inverse square of the uncertainty. The 
extrapolated zero-range energy is reported in the last 
row of column 7. The leading-order range dependence 
of Ezra ,2 (no) is cubic and we perform a fit of the form 
Co + C 3 ?~o + C 4 rg, weighted by the inverse square of the 
uncertainty. The extrapolated zero-range energy is re¬ 
ported in the last row of column 8. Table I shows 
that the zeroth-, first- and second-order extrapolation 
schemes yield zero-range energies that differ by at most 
0.0013E ho . This confirms that the range-dependence of 
the (3,2) ground state energy for the ro considered is well 
described by a Taylor series. Moreover, we note that 
Ezra, 2 (^ 0 ) for r 0 = 0.07a ho differs by only 0.0045E ho 
from the extrapolated zero-range energy. This suggests 
that Ezra, 2 (no) obtained at a single (relatively large) 
range provides a very good estimate for the zero-range 
energy. Circles, squares, and diamonds in Fig. 2 show 
the energies E ZRA ,o(no), E ZRA) i(r 0 ), and E ZRA , 2 (no), re¬ 
spectively. The fits (see the discussion above) are shown 
by lines. 

For systems with up to six particles (see the Supple¬ 
mental Material [27] for a summary), we believe that our 
basis sets for all ro are very close to complete. Specifi¬ 
cally, (i) the energy changes very little upon further en¬ 
largement of the basis set, (ii) the first- and second-order 
derivatives E^\ro) and E^ 2 ^(ro) are stable and their er¬ 
rorbars can be estimated reliably, (iii) the extrapolations 
of Ezra, o(jo), E Z RA,i(r 0 ), and Ezra^o) are in very 
good agreement, and (iv) the quantities 3E(r 0 ) — 4V[ r (r Q ) 
and Ezra, 1 (no) agree quite well [see the discussion af¬ 
ter Eq. (10)]. For systems with more than six particles, 
the construction of a nearly complete basis set is more 
challenging, especially for small r 0 . As an example, we 
discuss the (4,4) system; for the (4,3), (5,4), and (5,5) 
systems, the reader is referred to the Supplemental Ma¬ 
terial [27]. 



FIG. 3. (Color online) Ground state energy of the (4,4) sys¬ 
tem at unitarity. (a) Circles, squares, and diamonds show 
E Z ra,o()"o), E Z RA,i(r 0 ), and Ezra,2(^0), respectively, for the 
largest basis set considered. The errorbars of EzRA,o(ro) show 
the estimated basis set error AE(ro) (see column 3 of Ta¬ 
ble II); they extend below the data points but not above. The 
errorbars of Ezra, 1 (no) combine the estimated basis set error 
and the error of E^\ro) (see column 4 of Table II). Lastly, 
the errorbars of Ezra,2(^0) combine the estimated basis set 
error, and the errors of E^(ro) and E ( 2 ^(ro) (see column 5 
of Table II). Solid and dashed lines show the extrapolations 
of EzRA,o(ro) and EzRA,i(ro) to the zero-range limit, (b) 
Same quantities as in (a) but corrected for the estimated basis 
set errors. The open symbols show the energies Ezra, 0(no), 
Ezra,i()"o), an d Ezra,2(^0) reported in the last three columns 
of Table II. The uncertainty of AE(ro) is not accounted for 
by the errorbars. It can be seen that the basis set error low¬ 
ers the zero-range energy by about 0.05Eh o or, equivalently, 
0.4%. 


Table II summarizes our ECG results for the (4,4) sys¬ 
tem; the format is the same as that in Table I for the 
(3,2) system. The smallest range considered and the 
errorbars for the N = 8 system are larger than those 
for the N = 5 system. For a fixed ro, the quantities 
EW(r 0 ) and E^ 2 )(ro) for the N = 8 system are about 
twice as large as for the N = 5 system and their er- 
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TABLE II. Same as Table I but for the (4,4) system at unitarity. 


r o 

E ( r o) 

AE(r 0 ) 

E W (r 0 )^ 

^ho 

E (2) (r o)“t 

^ho 

^ZRA,0 ( r o) 

^zra,i ( r o ) e zra, 2 (H)) 

a ho 

E ha 

E hn 

E ha 

E ha E ho 

0.1 

12.329 

0.010 

2.07(10) 

-18(3) 

12.319 

12.113(10) 12.011(25) 

0.08 

12.287 

0.018 

2.44(16) 

-16(4) 

12.269 

12.073(13) 12.015(26) 

0.06 

12.230 

0.022 

2.72(25) 


12.208 

12.045(15) 

0.05 

12.204 

0.025 

2.71(32) 


12.179 

12.043(16) 

0.04 

12.184 

0.035 

2.56(55) 


12.149 

12.047(22) 

0 





12.015 

12.019 


rorbars are notably larger. The overall trends, however, 
are similar: (i) the energy E[tq) decreases with decreas¬ 
ing range, (ii) E^(ro) increases with decreasing range 
for ro > 0.06aho (for smaller ro, the trend reverses; we 
believe that this is a consequence of the numerics and 
not a real trend), and (iii) E^ 2 \r 0 ) becomes less nega¬ 
tive with decreasing range. Our numerics are not good 
enough to determine E^^ro) for r o < 0.06ah o - It can 
be seen, however, that the Ezra, 2 (A)) for ro = O.laho 
and 0.08aho agree quite well with the energies obtained 
by extrapolating -EzRA,o( r o) and E Z R A ,i( r o) to the zero- 
range limit (see the last row of Table II). Since the basis 
set error is non-negligible for N = 8, Fig. 3 shows the 
energies E ZRAfi (r 0 ) (circles), E ZRA ,i(ro) (squares), and 
-EzRA, 2 (ro) (diamonds) before correcting for the basis set 
error [Fig. 3(a)] and after correcting for the basis set error 
[Fig. 3(b)]. We extrapolate E ZRAt0 (r 0 ) and E ZRA s(r 0 ) 
for both the largest basis set considered [Fig. 3(a)] and 
the infinite basis set [Fig. 3(b)] to the zero-range limit by 
performing fits of the form cq + Cir 0 + C 2 ?’g and Co + C 2 rg, 
respectively (see solid and dashed lines in Fig. 3). To fit 
a function to -Ezra.i^o)) the data points are weighted 
by the inverse square of the uncertainty. For the (4,4) 
and larger systems, we do not fit to higher-order poly¬ 
nomials because (i) the number of data points is five or 
less and (ii) the errorbars are too large to determine the 
r[j dependence reliably. An alternative fit approach that 
includes the ?’g term is discussed in Appendix A. 

Table III summarizes our zero-range ground state en¬ 
ergies E(0) (column 2) obtained by extrapolating the en¬ 
ergies -Ezra.i^o): which have been shifted down by the 
estimated basis set error, to the zero-range limit. The 
errorbars given in parenthesis are estimated by combin¬ 
ing the zero-range extrapolation error, the uncertainty of 
the basis set error, and the uncertainty of E^^tq). 

Our ground state energy for the (2,1) system agrees ex¬ 
cellently with the semi-analytic energy obtained using the 
zero-range framework of Ref. [31]. For comparison, Ta¬ 
ble III includes the ground state energies from the litera¬ 
ture obtained by various methods. Column 3 of Table III 
reports the zero-range ground state energies Eecg calcu¬ 
lated by the ECG method in previous works [18, 20, 23]. 
Our results for the (2,2), (3,2) and (3,3) systems agree 
within errorbars with the literature results. For the (3,3) 
system, we provide a notably tighter errorbar. The ECG 
energy for the (4,4) system by Bradly et al. [23] is about 
0.9% higher than our (4,4) energy. Bradly et al. esti¬ 


mate that the error due to the use of a restricted ba¬ 
sis set is about 0.6% for the relative energy, translating 
to 0.53% for the total energy. This estimate is in rea¬ 
sonable agreement with the difference between their en¬ 
ergy and our energy. Column 5 reports the FN-DMC 
energies Edmci calculated for the square well potential 
with range O.Olciho [7]. The deviations between the FN- 
DMC and our ECG energies are due to the positive effec¬ 
tive range correction and the approximate nature of the 
nodal surface of the trial wave function (the latter domi¬ 
nates). Column 7 reports highly-improved FN-DMC en¬ 
ergies EdmC 2 [12]. These energies have been extrapo¬ 
lated to the zero-range limit. The FN-DMC energies 
from Ref. [12] agree very well with our ECG energies 
(the agreement is better than 0.6% and, for N = 8 and 
10, the errorbars overlap). Unfortunately, Ref. [12] con¬ 
sidered only spin-balanced systems. Columns 9 and 11 
report the energies EafmC 2 and EafmC 4 calculated us¬ 
ing the AFMC approach with q 2 and q 2 + q A dispersion 
relations, respectively [12], These energies have been ob¬ 
tained by applying a leading-order correction scheme to 
convert the finite lattice results to the infinite lattice limit 
but have not been extrapolated to the infinite lattice size 
limit [12, 32]. Note that the AFMC energies for fixed 
N but different dispersion relations do not agree within 
errorbars. The reason may be that the corrections due to 
the finite lattice spacing behave differently for the differ¬ 
ent dispersion relations and that the errorbars are purely 
statistical. Column 13 reports the configuration interac¬ 
tion (Cl) energies Eqi obtained using a limited Cl shell 
model space [13]. The authors of Ref. [13] noted that the 
two-body interaction strength was renormalized using an 
approach that could be improved upon. Improvement to 
both these aspects (enlarged Cl model space and refined 
renormalization approach) could change the Cl energies. 
Interestingly, the odd N Cl energies agree quite well with 
our ECG energies while the even N Cl energies are higher 
by between 1.3% and 3.2%. It is not clear to us what the 
origin of the different even and odd N behaviors is. Col¬ 
umn 15 reports the lattice MC energies lattice [10]. The 
lattice MC energies exhibit shell effects that are absent 
in the FN-DMC, AFMC, and—for small N —ECG ener¬ 
gies (our energies for N < 10 do not exhibit shell effects). 
The lattice MC energy for N = 8 is 3.2% lower than our 
ECG, reflecting the shell effects exhibited by the lattice 
MC energies in the small N regime. For systems with 
N > 10, the lattice MC energies are higher than or equal 
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TABLE III. Summary of our zero-range ground state energies at unitarity and comparison with literature results for systems 
with N < 10 and Ni — N a < 1. Column 2 reports the zero-range ground state energies E( 0) calculated in this work. Columns 
3 to 16 report ground state energies from the literature calculated by different methods and their percentage differences from 
E( 0). All energies are reported in units of E^ a . The ground state energy of the (2,1) system, obtained semi-analytically [31], 
is 4.272724i?ho. See text for more details. 


{N u N a ) 

m 

Eecg 

% Fx>MCi a % £r>MC2 b 

% 

-EaFMC2 C % 

-EaFMC4 C % 

E C i A 

% 

-^lattice 6 % 

(2,1) 

4.2727(1) 


4.281(4) 0.2 




4.279 

0.1 


(2,2) 

5.0091(4) 

5.0092(4) f 

5.051(9) 0.8 5.028(2) 

0.4 



5.138 

2.6 

5.071(+32/—75) 1.2 

(3,2) 

7.455(1) 

7.457(3) s 

7.61(1) 2.1 







(3,3) 

8.337(4) 

8.34(9) g 

8.64(3) 3.6 8.377(3) 

0.5 

8.26(1) -0.9 

8 .21(1) -1.6 

8.601 

3.2 

8.347(+80/—66) 0.1 

(4,3) 

11 .01(2) 


11.36(2) 3.2 




11.021 

0.1 


(4,4) 

12.02(3) 

12.13(l) h 

0.9 12.58(3) 4.7 12.04(1) 

0.2 

11.82(2) -1.7 

11.76(3) -2.2 

12.179 

1.3 

11.64(+11/-12) -3.2 

(5,4) 

15.24(9) 


15.69(1) 3.0 







(5,5) 

16.12(6) 


16.80(4) 4.2 16.10(1) 

-0.1 





16.05(+3/-7) -0.4 


a From Table II of Ref. [7]; the energies have been calculated for the square well potential with range ro = 0.01a; K1 , corresponding to 
r eff — 0.01a ho . 

b Read off from Fig. 4 of Ref. [12]; the FN-DMC energies have been extrapolated to the zero-range limit. 

c Read off from Fig. 4 of Ref. [12]; the errorbars only account for the statistical uncertainty. A leading-order correction scheme has been 
applied to convert the finite lattice results to the infinite lattice limit [12, 32], 
d From Table I of Ref. [13]; the Cl energies have been obtained for a finite shell model space and the two-body coupling constant has 
been renormalized by matching the two-particle ground state energy to the exact energy. 
e From Table VI of Ref. [10]; the upper and lower limits of the errorbars are different and separated by a slash. The errorbars account 
for statistical, fitting, finite volume, and spatial discretization errors, but do not account for systematic errors due to the contributions 
from excited states. We note that odd N systems were considered in Ref. [9]. The results in Ref. [9] were described as “preliminary” 
and are not included here. 

f From Ref. [20]; the energy has been obtained by solving the hyperangular Schrodinger equation. 

g From Table XXI of Ref. [18]; the energies have been extrapolated to the zero-range limit using the zeroth-order extrapolation scheme. 
h From Table II of Ref. [23]; the energy has been extrapolated to the zero-range limit using the zeroth-order extrapolation scheme. 


to (within errorbars) the FN-DMC energies £’dmC2 from 
Ref. [12]. The difference between the lattice MC ener¬ 
gies and FN-DMC energies for N > 10 is smallest for 
closed shell systems. Besides the results summarized in 
Table III, we also compared our ground state energies 
with DFT energies [14] for both even and odd N sys¬ 
tems. The DFT energies are 5% to 10% higher than our 
ECG energies. 

TABLE IV. Zero-range contact C(0) at unitarity for N = 3 — 
10. Column 2 reports the zero-range contact C(0) determined 
using the adiabatic energy relation. C(0) for the (1,1) system, 
obtained analytically from the implicit eigenequation derived 
in Ref. [33], is 4 a/2 na~° = 10.026513a“ o 1 . C(0) for the (2,1) 
system, obtained semi-analytically using the hyperspherical 
coordinate framework [31, 34, 35], is 10.468967a]]], 1 . 


(JVi, JV 2 ) 

C(0)a h o 

(2,1) 

10.469(1) 

(2,2) 

25.74(1) 

(3,2) 

25.20(1) 

(3,3) 

40.39(8) 

(3,4) 

38.2(2) 

(4,4) 

55.4(5) 

(5,4) 

56.9(9) 

(5,5) 

72.3(8) 


In addition to the energies, we calculate the contact at 
unitarity. To remove the leading-order range dependence, 
we analyze the quantities C'zra.o(fo) and Czra,i(fo)- 


While the energies -Ezra,o(fo) an< ^ -®zra,i(to) approach 
the ro = 0 limit from above and below, respectively, for 
all N considered, the contacts Czra,o(fo) and Czra,i(fo) 
approach the ?’ 0 = 0 limit from either above or below. 
Specifically, fitting Czra,o(fo) to a function of the form 
Co + Ciro + C 2 To, we find that Ci is positive for N = 4, 
very close to zero for N = 6, negative for N = 8, and 
again positive for N = 10. For the odd N systems, ci is 
always positive. The pair distribution functions exhibit 
an analogous range-dependence in the ro < r <C Oho re¬ 
gion (see Figs. 9 and 10 of the Supplemental Material for 
the N = 5 and 8 systems), suggesting that the intricate 
N and r 0 dependence of the contact is a real effect and 
not an artifact of our numerics. Our convergence studies 
support this interpretation. Table IV reports the zero- 
range contact C(0) for TV = 3 — 10 at unitarity obtained 
by extrapolating CzRA,i(ro) to the zero-range limit. The 
errorbars in parenthesis account for the zero-range ex¬ 
trapolation error and the basis set error. The ro = 0 ex¬ 
trapolations of C'zra,o(i’o) and of the contact extracted 
from the pair distribution functions agree with the val¬ 
ues reported in Table IV but have larger errorbars. The 
contact exhibits an interesting even-odd pattern. Specif¬ 
ically, for the N = 4 and 5 systems (and the 6 and 
7 systems, and the 8 and 9 systems), the contacts are 
roughly equal, reflecting the fact that these neighboring 
even-odd systems contain the same number of pairs. To 
zeroth-order, the contact scales as N2 times the contact 
of the two-body system, i.e., linearly with the number of 



pairs. Since (7(0) scales with N 2 , the r 0 <C r <C aho re¬ 
gion of the scaled pair distribution functions 4nPi 2 (r)r 2 
approximately collapse to a single curve if multiplied by 
N\. This approximate collapse is illustrated in Figs. 4(c) 
and 4(d). 



FIG. 4. (Color online) Panels (a) and (c) show the scaled pair 
distribution functions 47rPi2(r)r 2 and 4nNiPi2(r)r 2 , respec¬ 
tively, of the ground state at unitarity for the (2,2) system 
(solid line), (3,3) system (dashed line), (4,4) system (dotted 
line), and (5,5) system (dash-dotted line). Panels (b) and (d) 
show the scaled pair distribution functions 47rPi2(r)r 2 and 
4nNiPi2(r)r 2 , respectively, of the ground state at unitarity 
for the (2,1) system (solid line), (3,2) system (dashed line), 
(4,3) system (dotted line), and (5,4) system (dash-dotted 
line). The calculations are performed for ro = 0.06dho. 

Figure 5 shows the radial density P\(r) of the ground 
state at unitarity for even N and ro = 0.06ah o - We 
note that the convergence of the radial density in the 
small r regime is not as good as that of the pair distribu¬ 
tion function, especially for large N and small Vq. Pi{r) 
peaks at r = 0 for the (2,2) system, is relatively flat 
in the small r region for the (3,3) and (5,5) systems, and 
peaks around 0.6ah o for the (4,4) system. To estimate the 
range dependence, we calculate P\{r) for different ro for 
the (2,2), (3,3), and (4,4) systems. For a given system, 
the r < 0.5ah o region of P\{r) increases with decreas¬ 
ing two-body range r 0 (see Fig. 8 of the Supplemental 
Material [27]). The changes with ro are relatively small 
and the densities displayed in Fig. 5 show the generic 
behavior of trapped Fermi gases with short-range inter¬ 
actions. Figure 6 shows Pj(r), j = 1 and 2, for the odd 
N systems at unitarity for ?' 0 = 0.06ah o - Pi(r) and P 2 (r) 
peak at r = 0 for the (2,1) system, are relatively flat in 
the small r region for the (3,2) and (5,4) systems, and 
peak around 0.5ah o for the (4,3) system. We find that 
the range-dependence of the radial density for the odd 
N systems is similar to that for the even N systems (see 
Fig. 7 of the Supplemental Material [27]). 

To gain insights into the pairing of the particles, Fig. 7 
shows the integrated quantities Nj(r), 

Nj (r) = 4-irNj [ Pj^r^dr', (14) 

Jo 



FIG. 5. (Color online) Radial density Pi(r) of the ground 
state at unitarity for the (2,2) system (solid line), (3,3) system 
(dashed line), (4,4) system (dotted line), and (5,5) system 
(dash-dotted line). The calculations are performed for r o = 
0.06ah o . 



FIG. 6. (Color online) Panels (a) and (b) show the radial 
density of the majority species Pi(r) and the minority species 
P 2 (r), respectively, for the ground state of the (2,1) system 
(solid line), (3,2) system (dashed line), (4,3) system (dotted 
line), and (5,4) system (dash-dotted line). The calculations 
are performed for ro = 0.06aho- 


for the odd N systems. Solid and dashed lines show Nj(r) 
for the majority (j = 1) and minority (j = 2) species, 
respectively. Nj(r) monitors the number of particles of 
species j located between zero and r, and approaches 
Nj in the large r limit. We find that Ni(r) and N 2 (r) 
take, for N fixed, different values for all r, suggesting 
that there exists no core region where the systems are 
fully paired. This is in contrast to an earlier FN-DMC 
study [8], which suggested that the N = 9 system has a 
fully paired core. It should be noted that a fully paired 
core is expected in the large N limit [36]; however, how 
many particles are needed to be in the large N limit is 
not clear. 
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FIG. 7. (Color online) Solid and dashed lines show the in¬ 
tegrated quantities TVi(r) and N 2 (r), respectively, for odd N 
systems as a function of r. From bottom to top, the curves 
correspond to systems with N = 3, 5, 7, and 9. The hori¬ 
zontal dotted lines at 1 to 5 serve as a guide to the eye. The 
calculations are performed for ro = 0.06aho- 


IV. CONCLUSIONS 

This paper considered the ground state properties of 
trapped two-component Fermi gases at unitarity with 
up to ten particles. The calculations were performed 
for interspecies finite-range Gaussian interaction poten¬ 
tials using the ECG approach. Previous ECG calcula¬ 
tions were limited to TV = 3^6 and 8. The present 
work additionally considered the spin-imbalanced TV = 7 
and 9 systems with L n = 1“ symmetry and the spin- 
balanced TV = 10 system with L 71 ' = 0 + symmetry. A 
new range-correction scheme, which allows for the lead¬ 
ing and- in some cases—the sub-leading range depen¬ 
dence to be removed, was introduced. The accuracy of 
the range correction scheme was tested extensively for 
small TV systems (TV < 6) and then applied to larger sys¬ 
tems ( TV = 7—10). The resulting extrapolated zero-range 
energies have errorbars that range from 0.002% for TV = 3 
to 0.6% and 0.4% for TV = 9 and 10. The energies agree 
well with the FN-DMC energies from Ref. [12], suggest¬ 
ing that the zero-range energies of harmonically trapped 
two-component Fermi gases with TV < 10 (W — TV 2 = 0 or 
1) are now known with an accuracy better than 1%. The 
finite-range energies were reported for finite ro and all 
TV. These finite-range energies provide variational upper 
bounds and are expected to help assess the accuracy of 
future finite-range calculations (the range ro can be easily 
converted to the effective range). In addition to the en¬ 
ergy, the pair distribution functions and radial densities 
were analyzed. The Tan contacts obtained through the 
adiabatic and pair relations were found to agree within 
errorbars. 
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Appendix A: Additional comments on the 
range-correction scheme 

In the main text, we independently fit the quantities 
E Z ra,o(A)), -EzRA,i(ro)) and -Ezra, 2 (Vo)- The resulting 
zero-range energies were found to be in good agreement. 
This appendix discusses that a single correlated fit yields 
results that are consistent with those obtained from the 
independent fits. 

We assume that the ground state energy E(r 0 ) = 
E Z ra,o(u>) is a polynomial in the two-body interaction 
range r 0 , 

E(r 0 ) = c 0 + cir 0 + c 2 4 + c 3 4 + 0(4). (Al) 

Using Eq. (Al) to calculate EW(?'o) and E^ 2 \ro) and 
inserting the results into Eq. (10), we find 

-Ezra, i(T o) = c 0 - c 2 rl - 2c 3 4 + 0(4) (A2) 

and 

Ezra, 2 (^ 0 ) = c 0 + c 3 4 + 0(4)- (A3) 

As expected, the leading-order range-dependencies of 
Ezra,i(co) and Ezra, 2 (^ 0 ) are quadratic and cubic, re¬ 
spectively, and the functional forms of TTzrAj (j'o) are 
not independent. Specifically, the quadratic coefficient of 
Ezra,i(co) has the opposite sign but the same magnitude 
as that of E Z ra,o(A))j the cubic coefficient of Ezra, 1 (^ 0 ) 
has the opposite sign but twice the magnitude as that 
of EzRA,o( r o)j an d the cubic coefficient of Ezra, 2 (^ 0 ) is 
the same as that of E Z ra,o4o)- Interestingly, our in¬ 
dependent fits shown in Figs. 2 and 3 of the main text 
and Figs. 1-6 of the Supplemental Material are largely 
consistent with Eqs. (A1)-(A3). For example, our fits 
of Ezra,o4o) yield a negative 4 coefficient and those 
of Ezra, 1 (^ 0 ) yield a positive 4 coefficient. The magni¬ 
tudes of these coefficients, however, depend fairly sensi¬ 
tively on the number of terms included in the indepen¬ 
dent fits. 

As an alternative, we perform a simultaneous four- 
parameter fit of E ZRA ,o(r 0 ), Ezra.iM, and E z RA ,2 4o) 
using Eqs. (A1)-(A3). Each data point is weighted by the 
inverse square of the uncertainty [for Ezra, 0 (^ 0 ) we as- 
sume an uncertainty of 0.3AE(tq) and the uncertainties 
of Ezra, 1 (^ 0 ) an d E Z ra,24o) are given in Tables I and II 
of the main text and Tables I-VI of the Supplemental Ma¬ 
terial]. The resulting zero-range energies for TV = 3 — 10 
are 4.2726E ho , 5.0088£ ho , 7.454E ho , 8.335A ho , U.01A ho , 
12.027^0, 15.25-Eho 5 and 16.12 E ho: respectively. These 
energies lie within the errorbars of the zero-range ener¬ 
gies reported in Table III. The simultaneous fit yields a 
positive ci coefficient and negative c 2 and c 3 coefficients 
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for all TV. The C\ coefficient obtained from the indepen¬ 
dent fit of -Ezra,o(To) differs from that obtained from the 
simultaneous fit by less than 10% for all TV. 

We also apply the simultaneous fit approach to the 
contact. We fit our numerically obtained C'zRA,o( r o) and 


C’zra, 1(7*0) to functions of the form Co + CiTq + C 27 q + 
C 3 Tq and Co — C 2 Tq — 2c3Tq, respectively, for TV < 6, and 
to functions of the form co + Giro + C 2 r§ and cq — C 2 Tq, 
respectively, for TV = 7 — 10. The resulting zero-range 
contacts C(0) for TV = 3 — 10 lie within the errorbars of 
the zero-range contacts reported in Table IV. 
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TABLE I. Ground state energy of the (2,1) system at unitarity. Column 2 shows the finite-range energy for the largest basis 
set considered. The estimated basis set error AE(ro) is reported in column 3. Columns 4 and 5 report the quantities E^'ho) 
and E^(ro); errorbars are given in parenthesis. The energy derivatives are calculated for the largest basis set considered. 
Columns 6-8 report the energies Ezra, oho), EzRA.iho), and Ezra, 2 (^ 0 ). These energies account for the estimated basis set 
extrapolation error, i.e., E(ro) — AE(ro) is being used to calculate Ezra,.; ho) for j = 0, 1, and 2. The errorbars of Ezra, 1 ho) 
and Ezra, 2 ho) account for the uncertainties of E^\r 0 ) and E (2 \r 0 ) but do not account for the uncertainty of AE(ro). The 
last row reports the extrapolation of Ezra,; ho) to the zero-range limit. 


rg E(r 0 ) A E(r 0 ) p(l) / \ aho 

_ E ^_ ( ' 0> E ha 

0.12 4.29594 0.00001 0.1502(2) 

0.1 4.29276 0.00001 0.1683(5) 
0.08 4.28923 0.00001 0.1846(9) 
0.06 4.28538 0.00001 0.199(2) 

0.04 4.28128 0.00001 0.209(2) 

0.03 4.27916 0.00001 0.213(5) 

0.02 4.27702 0.00002 0.214(5) 

0.01 4.27488 0.00002 0.214(7) 

“0 


E (2) (r 0 ) 


-0.91(1) 

-0.87(1) 

-0.68(3) 

-0.58(3) 

-0.35(3) 

-0.23(5) 

-0.10(6) 

0.09(15) 


E ZRA.oTo) 

__ 

4.29593 

4.29275 

4.28922 

4.28537 

4.28127 

4.27915 

4.27700 

4.27486 

4.27266 


^ZHA.lIni) 
_ T 1., _ 

4.27791(2) 

4.27592(5) 

4.27445(7) 

4.27341(12) 

4.27290(8) 

4.27277(15) 

4.27271(10) 

4.27271(7) 

4.27271 


E ZRA,2(U)) 

_ Eh n _ 

4.27137(10) 

4.27155(10) 

4.27227(17) 

4.27237(17) 

4.27262(10) 

4.27266(17) 

4.27269(11) 

4.27272(8) 

4.27273 



FIG. 1. (Color online) Ground state energy of the (2,1) system at unitarity as a function of r 0 . Circles, squares, and diamonds 
show the energies Ezra, oho), EzRA.iho), and Ezra, 2 (^ 0 ), respectively, reported in the last three columns of Table I. The 
uncertainty of AE(ro) is not accounted for by the errorbars. Solid, dashed, and dotted lines show polynomial fits to Ezra, oho), 
EzRA,iho), and Ezra, 2 (^ 0 ). 
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TABLE II. Same as Table I but for the (2,2) system at unitarity. 


H) E ( r 0 ) 

Shn _Siu_ 

0.12 5.21941 
0.1 5.18697 
0.08 5.15329 
0.06 5.11846 
0.05 5.10068 
0.04 5.08271 
0.03 5.06455 
0.02 5.04623 
0 


A-EQo) 

Ehr> 


E w (r 0 ) 


a ho 

E br> 


0.00005 1.590(16) 

0.00006 1.654(12) 

0.00007 1.714(11) 
0.00010 1.766(16) 
0.00012 1.789(33) 

0.00013 1.808(22) 
0.00016 1.824(31) 
0.00026 1.834(31) 


£ (2) M 


-3.26(3) 

-3.03(8) 

-2.70(14) 

-2.22(15) 

-1.84(26) 

-1.37(32) 

-0.96(34) 

0.07(48) 


^ZRA.of^o) 

_ E h' l _ 

5.21936 

5.18691 

5.15322 

5.11836 

5.10056 

5.08258 

5.06439 

5.04597 

5.00874 


-EzRA.l(m) 

_ ^'ll _ 

5.02861(192) 

5.02149(120) 

5.01611(88) 

5.01238(96) 

5.01113(165) 

5.01027(88) 

5.00969(93) 

5.00930(62) 

5.00912 


-Ezra, 2 C 0 ) 

_ Eh.. _ 

5.00515(214) 

5.00637(160) 

5.00746(133) 

5.00838(123) 

5.00883(198) 

5.00918(114) 

5.00926(108) 

5.00931(72) 

5.00939 



FIG. 2. (Color online) Same as Fig. 1 but for the (2,2) system at unitarity. 


TABLE III. Same as Table I but for the (3,3) system at unitarity. 


rp 

Ahn 

E(j- 0 ) 

Ehn 

AE(r 0 ) 

Ehn 

°l 2 

<3 |K) 

O 

K] 

e (2) Mlfc 

^ZRA,0 ( r o) 
E h „ 

^ZRA.I^O) 

Ehn 

^ZRA,2 (^o) 

Ehn 

0.1 

8.6097 

0.0009 

2.249(8) 

-11(1) 

8.6088 

8.3839(8) 

8.3289(58) 

0.08 8.5628 

0.0012 

2.491(14) 

-10(1) 

8.5616 

8.3623(11) 

8.3303(43) 

0.07 8.5380 

0.0018 

2.597(16) 

-9(1) 

8.5362 

8.3544(11) 

8.3324(36) 

0.06 

8.5105 

0.0021 

2.697(21) 

-8(2) 

8.5084 

8.3466(13) 

8.3322(49) 

0.05 

8.4839 

0.0023 

2.779(35) 

-6(2) 

8.4816 

8.3427(18) 

8.3352(43) 

0.04 

8.4558 

0.0028 

2.789(62) 


8.4530 

8.3415(25) 


0.03 8.4281 

0.0034 

2.791(110) 


8.4247 

8.3410(33) 


0 





8.3413 

8.3369 

8.3374 


TABLE IV. Same as Table I but for the (4,3) system at unitarity. 


rp 

a ho 

E(r 0 ) 

Ehn 

AE(r 0 ) 

E ha 

E ll) (r 0 )^ £ ( 2 ) Mf^ 

^ZRA,0 ( r o) 

^ho 

^ZRA, l(^o) ^ZRA,2(^o) 

0.1 

11.166 

0.007 

0.94(3) -11(2) 

11.160 

11.065(3) 11.010(13) 

0.08 

11.146 

0.009 

1.18(12) -9(2) 

11.137 

11.043(10) 11.014(16) 

0.06 

11.120 

0.011 

1.33(28) 

11.108 

11.029(17) 

0.04 

11.095 

0.017 

1.13(50) 

11.078 

11.033(20) 

0 




11.004 

11.014 
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FIG. 3. (Color online) Same as Fig. 1 but for the (3,3) system at unitarity. 



FIG. 4 . (Color online) Ground state energy of the ( 4 , 3 ) system at unitarity. (a) Circles, squares, and diamonds show Ezra.o^o), 
F?ZRA,i(ro), and EzRA,2(ro), respectively, for the largest basis set considered. The errorbars of Ezra.o^o) show the estimated 
basis set error AE(ro) (see column 3 of Table IV); they extend below the data points but not above. The errorbars of Ezra,i (to) 
combine the estimated basis set error and the error of E ( 1 ) (ro) (see column 4 of Table IV). Lastly, the errorbars of Ezra,2(^0) 
combine the estimated basis set error, and the errors of E^(ro) and E^(r 0 ) (see column 5 of Table IV). Solid and dashed 
lines show the extrapolations of EzRA,o(fo) and EzRA,i(ro) to the zero-range limit, (b) Same quantities as in (a) but corrected 
for the estimated basis set errors. The open symbols show the energies Ezra,o(B)), EzRA,i(ro), and Ezra,2(^0) reported in the 
last three columns of Table IV. The uncertainty of AE(ro) is not accounted for by the errorbars. 


TABLE V. Same as Table I but for the (5,4) system at unitarity. 


7*0 

E (r 0 ) 

AE(r 0 ) 

EW(ro)fb 

E (2) (r 0 )^ 

^ZRA,0 ( r o) 

B ZRA,l( r o) -EzRA,2( r o) 

a ho 

E ho 

E ho 

E ho 


0.12 

15.622 

0.032 

2.09(17) 

-16(3) 

15.590 

15.340(20) 15.224(42) 

0.1 

15.584 

0.045 

2.37(18) 

-13(4) 

15.539 

15.302(18) 15.237(38) 

0.08 

15.536 

0.056 

2.59(31) 


15.480 

15.273(25) 

0.06 

15.490 

0.067 

2.14(85) 


15.423 

15.295(51) 

0 





15.227 

15.236 
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FIG. 5. (Color online) Same as Fig. 4 but for the (5,4) system at unitarity. 


TABLE VI. Same as Table I but for the (5,5) system at unitarity. 


r 0 

a ho 

E(r 0 ) 

E ha 

A E(r 0 ) 

E ho 

°l 2 

vk 

o 

<N 

K) 

A 
<3 |K] 

o 

£ ZRA,o( r o) 

E hn 

^ZRA,1 ( r 0 ) E ZRA,2( r o) 
E ha E ho 

0.12 

16.659 

0.024 

2.64(14) —28(6) 

16.635 

16.318(17) 16.127(60) 

0.1 

16.594 

0.030 

3.17(25) -24(4) 

16.564 

16.247(25) 16.117(45) 

0.08 

16.522 

0.037 

3.55(32) 

16.485 

16.201(26) 

0.06 

16.439 

0.040 

3.40(65) 

16.399 

16.195(39) 

0 




16.095 

16.124 



FIG. 6. (Color online) Same as Fig. 4 but for the (5,5) system at unitarity. 


























5 



FIG. 7. (Color online) Panels (a) and (b) show the radial densities Pi(r) and /^(r) of the majority and minority species, 
respectively, for the ground state of the (3,2) system at unitarity for ro = 0.06aho (dotted line), 0.04ah o (dashed line), and 
0.02ah o (solid line). The insets of panels (a) and (b) show blow-ups of the small r regions. 



FIG. 8. (Color online) Radial density Pi(r) for the ground state of the (4,4) system at unitarity for ro = O.laho (dotted line), 
0.08ah o (dashed line), and 0.06ah o (solid line). The inset shows a blow-up of the small r region. 



FIG. 9. (Color online) Panel (a) shows the scaled pair distribution function 47rPi2(r)r 2 of the ground state for the (3,2) system 
at unitarity for ro = 0.06ah o (dotted line), O.tMdho (dashed line), and 0.02ah o (solid line). Panel (b) shows a blow-up of the 
small r region. 
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r/a ho 

FIG. 10. (Color online) Panel (a) shows the scaled pair distribution function 4jrPi2(r)r 2 of the ground state for the (4,4) 
system at unitarity for r o = O.laho (dotted line), 0.08ah o (dashed line), and 0.06ah o (solid line). Panel (b) shows a blow-up of 
the small r region. 





